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SUMMARY 


The steady-state incompressible Navier-Stokes equations in two dimensions are solved numer- 
ically using the artificial compressibility formulation. The convective terms are upwind-differenced 
using a flux-difference split approach that has uniformly high accuracy throughout the interior grid 
points. The viscous fluxes are differenced using second-order-accurate central differences. The 
numerical system of equations is solved using an implicit line-relaxation scheme. Although the 
current study is limited to steady-state problems, it is shown that this entire formulation can be 
used for solving unsteady problems. Characteristic boundary conditions are formulated and used 
in the solution procedure. The overall scheme is capable of being run at extremely large pseudotime 
steps, leading to fast convergence. Three test cases are presented to demonstrate the accuracy and 
robustness of the code. These are the flow in a square-driven cavity, flow over a backward-facing 
step, and flow around a two-dimensional circular cylinder. 

INTRODUCTION 


The motivation for the current work is a desire to find an efficient method of solution for the 
incompressible Navier-Stokes equations for complex three-dimensional (3-D) geometries. Several 
methods exist which reformulate the equations into a nonprimitive-variable form. These include 
a velocity-vorticity method as proposed by Fasel [1], and a vector potential-vorticity method pro- 
posed by Aziz and Heliums [2]. Each of these methods requires the solution of three Poisson equa- 
tions at each time level, which is one reason why these methods have not become popular. Further 
work on the velocity-vorticity method and the use of a direct solver is being done by Hafez, Da- 
cles, and Soliman [3]. These methods currently appear to be too expensive for solving large 3-D 
problems. This leaves methods formulated in primitive variables. The primary difficulty in solving 
the incompressible Navier-Stokes equations in primitive variables stems from the lack of a time 
derivative in the continuity equation. There is no straightforward way to iteratively march these 
equations in time and ensure a divergence-free velocity field. The marker-and-cell (MAC) method 
of Harlow and Welch [4] accomplishes this by solving for the velocity field from the momentum 
equations and for the pressure field from a Poisson equation which is derived by taking the diver- 
gence of the momentum equations. This method can be very costly because of the need to solve 
the Poisson equation at each time iteration, and because there is no strong coupling between the 
velocity field and the pressure field. 

Another approach which has been used extensively is that of the fractional-step method origi- 
nally introduced by Chorin [5], and used by Kim and Moin [6], and Rosenfeld, Kwak, and Vinokur 
[7]. This method advances the solution in time using two (or more) steps. First the momentum 
equations are solved for an intermediate velocity field which will generally not be divergence-free. 
In the next step, a pressure field is obtained which will map the intermediate velocity field into a 
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divergence-free field, thus obtaining the solution at the next time level. The second step generally 
requires the solution of a Poisson equation in pressure. 

The method chosen for the current work is known as the pseudocompressibility method and 
was first introduced by Chorin [8]. It has been used with much success by Kwak et al. [9] for 
solving complex incompressible flow problems in generalized coordinates. In this formulation, 
a time derivative of pressure is added to the continuity equation. Together with the momentum 
equations, this forms a hyperbolic system of equations which can be marched in pseudotime to a 
steady-state solution. The method can also be extended to solve time-dependent problems [10,11] 
by using subiterations in pseudotime at each physical time step. If all that is desired is the steady- 
state solution to a problem, the pseudocompressibility method can be a very efficient formulation 
because it does not require that a divergence-free velocity field be obtained at each iteration. The 
addition of the time derivative of pressure to the continuity creates a hyperbolic system of equations 
complete with artificial pressure waves of finite speed. When the solution converges to a steady 
state, a divergence-free solution is obtained. Since this is the case, many of the well-developed 
compressible flow algorithms can be used for this method. 

Many previous applications of this method have used central differencing of the convective 
fluxes. This approach also requires that artificial dissipation be explicitly added in order to damp 
out the spurious oscillations which are a result of the nonlinearity of the convective fluxes. Such a 
scheme can be difficult to apply because the artificial dissipation parameters must be adjusted for 
each specific calculation. The use of the artificial dissipation will also tend to hamper the accuracy 
of the calculations [12]. To avoid the problems associated with central differencing, an upwind- 
differencing scheme is considered here. Of most recent interest has been the class of upwind- 
differencing schemes which bias the differencing based on the sign of the eigenvalues of the con- 
vective flux Jacobians. A number of these types of schemes have been developed in conjunction 
with solving the Euler equations and the compressible Navier-Stokes equations [13-16]. The de- 
velopment of these schemes base the upwind differencing on the physics of the Riemann problem. 
In the case of the pseudocompressibility method, the upwind differencing is merely a way of using 
the physics of the artificial waves to obtain a smooth numerical solution. 

Much of the current development of upwind-differencing schemes has focused on the ability to 
resolve sharp discontinuities without kinks or overshoots. By limiting the order of the differencing 
at points near the discontinuities, and thereby increasing the dissipation provided by the differ- 
encing, the schemes have the total variation diminishing (TVD) property. Applications of TVD 
schemes to the incompressible equations were done by Hartwich and Hsu [17,18], and Gorsfci [19]. 
These investigators were able to obtain 3- D solutions which were third-order-accurate in the con- 
vective terms, except near regions of large change in gradient, where the order of the differencing 
was reduced to increase the amount of dissipation added. 

Since solutions to the incompressible equations do not have strong discontinuities such as 
shocks, it is reasoned that the incompressible equations could be solved without the need for any 
limiting, and that flux-difference splitting of uniformly high order could be used. This paper at- 
tempts to show that this is so by using a flux-difference splitting type of formulation similar to that 
used for compressible flow in [13,14]. The current work concentrates on developing this scheme 
with the use of a two-dimensional (2-D) flow solver using fifth-order upwind differencing of the 
convective terms. Since the development of the upwind-differencing schemes considered here is 
based upon an analysis of a one-dimensional (1-D) hyperbolic conservation law, the use of a 2-D 
code for the initial testing done here will not be out of line from the desired goal of a 3-D al- 
gorithm. This will expedite much of the code development because of the smaller computational 
requirements of a 2-D code and because of the relative ease with which 2-D results can be analyzed, 
compared, and presented. 
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In the following sections, the details of the 2-D code are presented, including the govern- 
ing equations and the similarity transformation for the Jacobian matrix of the convective fluxes. 
The specific details of the upwind scheme are given, followed by the details of the implicit line- 
relaxation scheme used to solve the equations. Some boundary conditions based on the method of 
characteristics have been developed, and are presented. The computed results section shows the 
robustness and accuracy of the code by presenting three sample problems — the flow inside a driven 
cavity, the flow over a backward-facing step, and the flow over a circular cylinder. 

This work was partially sponsored by NASA Marshall Space Flight Center. 


GOVERNING EQUATIONS 


The equations governing constant-density viscous flow are presented here in nondimensional 
form. Following the pseudocompressibility formulation, a time derivative of pressure is added to 
the continuity equation resulting in 



/ du dv \ 
\3x + dy) 


= 0 


( 1 ) 


where u and v are velocity components in the x and y directions, respectively, and jS is known as 
the pseudocompressibility constant. Here, r represents pseudotime and is not related in any way 
to physical time. Combining equation (1) with the momentum equations for the incompressible 
Navier-Stokes equations results in the following system in Cartesian coordinates 
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where p is the pressure. In this formulation the Reynolds stress has been approximated as a function 
of the strain-rate tensor, and thus v represents a sum of the kinematic viscosity and the turbulent 
eddy viscosity. The equations in (2) are transformed into generalized curvilinear coordinates given 
by 

t = Z(x,y) 

V = ri(x,y) 

The equations are then given by 
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where J is the Jacobian of the transformation, the metrics of the transformation are 


. dr\ 


and the convective fluxes arc given by 
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where the contravariant velocity components, U and V are defined as 
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In deriving the viscous fluxes, constant viscosity was assumed for simplicity and because 
initially only laminar flow calculations are being performed. This simplification is not necessary 
and will be removed in the future. The viscous fluxes are given by 
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where v is the kinematic viscosity. If an orthogonal grid is assumed, the viscous fluxes reduce to 
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The upwind scheme requires the use of the eigensystem of the Jacobian matrix of the con- 
vective flux vectors. The 2-D eigensystem is presented here. For the 3-D equations, see Rogers, 
Chang, and Kwak [20], or Hartwich and Hsu [17]. Beware, however, that the transformation given 
by the latter can become singular for certain values of metrics. 


The generalized flux vector for the 2-D system of equations is given by 

A 
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where = E,F for i = 1, 2, respectively, and the metrics are represented by 
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where f i = = V> and the scaled contravariant velocity component is 


Q = k x u + k y v 

The Jacobian matrix of this generalized flux vector is given by 

fiky 
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A similarity transform for the Jacobian matrix is derived here of the form 

Ai = XiAiX- 1 

where 

-Q 0 0 ' 

A, = 0 Q + c 0 

.0 0 Q — c. 

and where c is the scaled artificial speed of sound given by 
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This variable will always be greater than Q, so the second eigenvalue will always be positive and 
the third eigenvalue will always be negative. The matrix of the right eigenvectors is given by 
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and its inverse is given by 
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UPWIND DIFFERENCING 


Upwind differencing is used to numerically compute the convective-flux derivatives. The 
upwind-differencing scheme is derived from 1-D considerations, and then is applied to each coor- 
dinate direction separately. Flux-difference splitting is used here to bias the differencing based on 
the sign of the eigenvalues of the convective-flux Jacobian. The scheme as presented here was orig- 
inally derived by Roe [14] as an approximate Riemann solver for the compressible gasdynamics 
equations. 

The derivative of the convective flux in the £ direction is approximated by 



where Ej+ 1/2 is a numerical flux and j is the discrete spatial index for the £ direction. 

The numerical flux is given by 

E j+ 1/2 = \ [E(Dj+i) + E(Dj) - <£ ;+ i /2 ] ( 14) 

For 0 ;+ , ji =0 this represents a second-order, central-difference scheme. The <A ;+1 / 2 is a dissipation 
term. A first-order upwind scheme is given by 


<£;'+i/2 - AEf +i /2 — AEj +1 / 2 (15) 

where A E ± is the flux difference across positive or negative traveling waves. The flux difference 
is computed as 

A^ 1/2 = ^(i))AD ; , 1/2 (16) 

where the A operator is given by 


A Dj+\/2 = Dj+ 1 - Dj 

The plus (minus) Jacobian matrix has only positive (negative) eigenvalues and is computed from 

A ± = X x A±X x x 

1 (17) 

A* = ^(Ai ± |Aj |) 

This Jacobian matrix is evaluated using some intermediate value which is a function of the sur- 
rounding points, j and j + 1. The Roe properties [14], which are necessary for a conservative 
scheme, are satisfied if this is taken as the average of the surrounding values. Thus 

D=j(D j+1 + Dj) (18) 


A scheme of arbitrary order may be derived using these flux differences. Implementation 
of higher-order-accurate approximations in an explicit scheme do not require significandy more 
computational time if the flux differences A E ± are all computed at once for a single line. A third- 
order upwind flux is defined by 

<£7+1/2 = - j(A£)/_i/2 - A£ ; + +1 /2 + A£ ;+1 /2 - A£ ;+3 / 2 ) ( 19 ) 


The primary problem with using schemes of accuracy greater than third order occurs at the 
boundaries. Large stencils will require special treatment at the boundaries, and a reduction of order 
is necessary. Therefore, when going to a higher-order-accurate scheme, compactness is desirable. 
Such a scheme was derived by Rai [21] using a fifth-order-accurate, upwind-biased stencil. A fifth- 
order, fully upwind difference would require 1 1 points, but this upwind-biased scheme requires only 
seven points. It is given by 

<£/+i/2 = ~3q[ ~ 2A Ej_ 3 /2 + 1 1 A Ej_i/2 ~ 6&Ej +l / 2 ~ 3A£' ; + +3/2 ^ 20 ) 

+ 2AE }+5 /2 — HAEj+3/2 + 6Ai? ;+1 /2 + 3A£7 _j^ 2 ] 
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Next to the boundary, near- second-order accuracy can be maintained by the third- and fifth- 
order schemes by using the following 


<t>j+\/2 - * 


AE j+i /2 ~ AE j+ 1/2 


( 21 ) 


For e - 0, this flux leads to a second-order central difference. For e = 1, this is the same as the 
first-order dissipation term given by equation (15). By including a nonzero e, dissipation is added 
to the second-order, central-difference scheme to help suppress any oscillations. A value of e = 
0.01 is used for all of the results presented in this paper. 


IMPLICIT SCHEME 


This section describes the way in which equation (3) is numerically represented and solved. 
Application of a first-order backward Euler formula to this system of equations yields the following 
delta-form equation 
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/At 


1 + 



(D ** 1 - D n ) = -R n 
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( 22 ) 


where the superscript n is the pseudotime iteration count and the vector R is the residual vector. 
When the formula given in equation (13) is followed and a second-order, central-difference formula 
is applied to the viscous terms, at a point Xj <k , Vj,k the numerical approximation to the residual vector 
is given by 


Rj,k 
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A( A t) 

(•fi?v).7+l ,k (^ti)/— l,fc ( Fv)j,k+l ( Fy)j,k—l 


(23) 
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where the numerical fluxes E and F are evaluated using equation (14) with either the first-, third-, or 
fifth-order dissipation term given in equations (15, 19,20), respectively. The generalized coordinates 
are chosen so that A £ and A 77 are equal to one. 


The formation of the exact Jacobian matrix of the residual vector will be too expensive for 
practical consideration, particularly when higher-order upwind differencing is used, so the implicit 
side formulation will be limited to using the residual resulting from the first-order upwind differ- 
encing. By applying the first-order dissipation term in equation (15) to the convective terms, the 
residual is given by 
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The exact Jacobian matrix of the residual vector will form a banded matrix of the form 


dR 

dD 


= B 


« QRjJi dRj^k dRjj± r) ft dRj,k 

[dD j>k ^ dDj- hk ’ dD,y dD j+hk ’ V ’ dD jM ! 


(25) 


7 



where B refers to a banded matrix. By using approximate Jacobians of the flux differences as 
derived and analyzed by Barth [22], the implicit side of the numerical equation is formed using the 
following terms 
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where A = Ai , B - Ai in equation (8), and 

A* = XiAfXf 1 
= X 2 A±X? 

Only the orthogonal mesh terms are retained for the implicit viscous terms resulting in 
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The matrix I m is a modified identity matrix given by 


(26) 
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The numerical system of equations thus formed is solved using a line-relaxation method. In 
this procedure, the entire numerical matrix equation is first formed from values at the previous time 
level. At this point the numerical equation is stored as a banded matrix of the form 

B[V,0 f ...,0,X,Y,Z,0,- >0,W]AD = R 

where A D = D n+1 — D n and V,W,X, Y, and Z are vectors of 3 by 3 blocks which lie on the diag- 
onals of the banded matrix, with the Y vector on the main diagonal. This matrix is approximately 
solved using an iterative approach. One family of lines is used as the sweep direction. Using, 
for example, the £ family, a tridiagonal matrix is formed by multiplying the elements outside the 
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tridiagonal band by the current A D and shifting them over to the right-hand side. This can be 
represented by the following 

BIX,Y,Z](AD) 1 " = R- AD‘ tik _,V - AD l jMt W 

where l is an iteration index. This system can be solved most efficiently by first performing and 
storing the LU decomposition of the tridiagonal matrix before the iteration is begun. Then for 
each iteration, the right-hand side is formed using the latest known A D, and the entire system is 
back-solved. The LU decomposition can be entirely vectorized, but the backsolution is inherently 
recursive and cannot be vectorized. 


UNSTEADY FORMULATION 


The current scheme can be easily extended to solve unsteady problems with the use of subiter- 
ations in pseudotime at each physical time step. The details of this formulation are given by Rogers 
and Kwak [10] and by Athavale and Merkle [11], and a short summary of this follows. First, the 
time derivative in the momentum equations is discretized using a second-order backward Euler 
formula, resulting in 

— (1.5 D^ 1 -2Z> n + 0 .5D*- 1 ) = -IT* 1 
At v ' 

where R is the same residual vector as in equation (23). Here physical time is denoted by t and 
the superscript n denotes the solution at time t = nA t. This equation leaves no way to update the 
pressure to the next time level because of the I m matrix on the left-hand side. Here the continuity 
equation is replaced with an artificial compressibility relation and a pseudotime level is introduced, 
resulting in 


I tT _ J D n+1 - m ) = -2D n + 0.5D n - 1 ) 

ZA Z 


Here pseudotime is denoted by t and the superscript m represents a subiteration index in pseudo- 
time. The matrix I tT is a diagonal matrix given by 


Itr = diag 


J_ hi LL 

At’ At ’ At 


After linearizing the residual about the n+ 1 , m time level, the following equation is obtained 



(28) 


It can be seen that this equation is very similar to its steady-state counterpart. The additional right- 
hand side terms and the different diagonal matrix on the left-hand side are the only differences 
between the two. This makes it quite simple to program a code capable of using this scheme to 
solve both unsteady and steady-state problems. 
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BOUNDARY CONDITIONS 


Implicit boundary conditions are used at all of the boundaries; this helps make possible the 
use of large time steps. At a no- slip surface, the velocity is specified to be zero, and the pressure at 
the boundary is obtained by specifying that the pressure gradient normal to the wall be zero. The 
boundary conditions used for inflow and outflow regions are based on the method of characteristics. 
The formulation of these boundary conditions is similar to that given by Merkle and Tsai [23], but 
the implementation is slightly different. The scheme is derived here for a £ = constant boundary, 
with similar results for a rj = constant boundary. The finite-speed waves which arise with the use 
of artificial compressibility are governed by the following 

dD _ dE _ dE dD _ , dD _ j dD 

dr dZ dD dZ dZ d( 


then 



= -AX" 1 


dD 

dZ 


(29) 


If one were to move the X -1 matrix inside the spatial and time derivative, then it can be seen that 
this would be a system of scaler equations, each with the form of a wave equation. The sigjt of 
the eigenvalues in the A matrix determines the direction of travel of the wave. For each positive 
(negative) eigenvalue, there is a wave propagating information in the positive (negative) £ direction. 
The number of positive or negative eigenvalues determines the number of characteristic waves 
propagating information from the interior of the computational domain to the boundary. Thus, at 
the boundary we will use these characteristics which bring information from the interior as part of 
our boundary conditions. The rest of the information should come from outside the computational 
domain, and we are free to specify some variables. 


There will be either one or two characteristics traveling toward the boundary from the interior 
because there is always at least one positive eigenvalue and one negative eigenvalue. In order to 
select the proper characteristic waves, equation (29) is multiplied by a diagonal selection matrix L 
which has an entry of 1 in the position of the eigenvalue we wish to select, and zeros elsewhere. 
Thus 


LX" 1 


dD 

dr 


= -LAX" 1 


dD 

dZ 


(30) 


Replacing the time derivative with an implicit Euler time step gives 


/ LX" 1 
\ /At 


+ iAX ~'|) 


(r 1 


— D n ) = -LAX" 1 


dD n 

dZ 


(31) 


This gives either one or two relations, depending on the number of nonzero elements in L. To 
complete the set of equations, some variables must be specified to be constant. Here is defined a 
vector Q of the variables to be held constant, such that 


d Q 
dr 


= 0 


dn dD 
dD dr 


dn 

dD 


= 0 


Combining equations (31) and (32) gives 


/ LX" 1 
\ /At 
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(32) 


(33) 
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Equation (33) can be used to update the variables implicitly at any of the inflow or outflow bound- 
aries with the proper choice of L and £2 . 


Inflow Boundary 


At the inflow, there will be one characteristic wave traveling out of the computational domain 
since fluid is traveling into the domain. If the incoming fluid is traveling in the positive £ direction, 
then 

Q > 0 

Q + c > 0 


Q - c < 0 

This third eigenvalue will be the one we wish to select, and so L will have a 1 for the third diagonal 
entry. If the incoming fluid is traveling in the negative £ direction, then 


Q < 0 
Q + c > 0 
Q - c < 0 


and the second eigenvalue is the one corresponding to the wave propagation out of the computa- 
tional domain, requiring a 1 in the second diagonal entry of L. 


Two different sets of specified variables have been used successfully for inflow boundaries. 
One set consists of the total pressure and the cross-flow velocity. This set is useful for problems in 
which the inflow velocity profile is not known. For this set the Q vector is 


-p+ j(u 2 + v 2 ) ■ 

dCl 

-1 

u 

v ' 
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0 

0 

0 

V 

’ dD ~ 

.0 

0 

1 . 


The second possible set of specified variables consists of the velocity components. This is useful 
for problems in which a specific velocity profile is desired at the inflow. The Cl vector for this is 


-0- 

da _ 

■o 

0 

o- 
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0 
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. V - 

; dD ~ 

.0 
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1 . 


Outflow Boundary 


At the outflow boundary there are two characteristic waves traveling out of the computational 
domain since fluid is also leaving the domain. If the fluid is traveling along the positive £ direction, 
then 

Q> 0 
Q + c> 0 
Q - c < 0 


11 



and we require a 1 in the first two diagonal entries of the L matrix. If the fluid is traveling in the 
negative £ direction, then 

Q< 0 
Q+c>0 

Q - c <0 


and we require a 1 in both the first and third diagonal entries of the L matrix. 


For all of the test problems presented in this paper, static pressure was specified at the outflow 
boundary, resulting in 


■p' 


■1 

0 
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.0. 

; ~dD ~ 
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0 

0. 


COMPUTED RESULTS 


The code was run for three 2-D laminar-flow test cases. These are a driven square-cavity 
flow, flow over a backward-facing step, and flow over a circular cylinder. These cases were chosen 
because they each have been studied previously by others in experiments or in computations. The 
computing times reported here are the CPU seconds used on a Cray 2. For comparison, these 
times are nearly the same as those obtained running on a Cray XMP-48. The computations are run 
until the maximum residual has converged over six orders of magnitude, the maximum divergence 
of velocity over all the points is less than 10 " 4 , and the flow quantities being measured have 
approached a steady-state value in at least four significant digits. 

For each of the test cases presented, the larger the time step A r, the better the convergence 
was, provided the solution remained stable. In all of the cases presented here the solution remained 
stable no matter how large a time step was used, so the time step was set to 10 12 , which effectively 
reduced the 1 / A r term to zero. The choice of /3 for each case was arrived at through numerical 
convergence tests. It was found that the convergence was quite sensitive to the value of / 3 , and in 
some cases, the choice of /3 could cause the scheme to become unstable. For most cases, however, 
it was easy to find a range of /? for which the code would converge very quickly. The convergence 
of the current formulation is degraded by the errors introduced by the approximate Jacobians on 
the left-hand side of the equations and by the fact that the whole system of equations is not exactly 
solved by the line-relaxation process. If it were possible to use the exact Jacobians and solve the 
system exactly, then this would be a Newton iteration, in which case one would expect to nave 
quadratic convergence when using a very large time step for any value of f3. Analysis of these 
errors and their relationship to (3 is under way, and it is hoped that a guideline for choosing /3 and 
for minimizing the eigenvalues of the amplification matrix can be obtained. Until such a guideline 
is found, the numerical tests will have to suffice. 


Driven Cavity Flow 


The 2-D flow in a driven square cavity whose top wall moves with a uniform velocity has been 
used rather extensively as a validation test case by several authors in the recent past. It provides 
a good test case in that there is no primary flow direction and the boundary conditions are very 


12 



simple to use. Ghia, Ghia, and Shin [24] presented extensive numerical data obtained from their 
multigrid vorticity-stream function formulation using very fine grids. They reported results which 
agreed well with other computational efforts. Other recent computational work involving this par- 
ticular geometry include Schreiber and Keller [25] who use a vorticity-stream function formulation; 
Kim and Moin [6] who use a fractional-step method in primitive variables in conjunction with ap- 
proximate factorization; Vanka [26] who uses a multigrid technique in primitive variables; and 
Benjamin and Denny [27] who use a centrally differenced vorticity-stream function formulation in 
conjunction with an ADI scheme. 

The current calculations attempt to maintain the accuracy of these authors while using fewer 
grid points. The flow is calculated for Reynolds numbers of 100, 400, 1000, 3200, 5000, 7500, and 
10,000 using a grid of 81 by 81 points where the points are clustered toward the walls. This grid 
is shown in figure 1. The value of the artificial compressibility /? was set to 20 for the Reynolds 
number of 100, to 10 for the 400 Reynolds number ease, to 2 for a Reynolds number of 1000, and to 
1 for the higher Reynolds numbers. The implicit line relaxation used 11 sweeps in the f -direction 
for each iteration. 

The velocity components on the lines passing through the geometric center of the cavity are 
compared to the results of Ghia, Ghia, and Shin [24] in figure 2. The u- velocity component is plotted 
along the y-axis for the different Reynolds numbers in figure 2a. The origins of the plots have been 
shifted to the left for each successive Reynolds number case. The data of Ghia were obtained from 
a uniform grid of 129 by 129 points for Reynolds numbers up to 3200, and a uniform grid of 257 
by 257 points for Reynolds numbers 5000 and greater. It is noted that these two computed results 
agree very well with each other. In figure 2b, the v-velocity component is plotted along the x-axis 
passing through the geometric center for the different Reynolds numbers, llie origins of these plots 
are shifted up for each successive Reynolds number case. Again, good agreement is seen between 
the two computed results. 

In table 1, the stream function and vorticity quantities are given for the core of the primary 
vortex for all the Reynolds numbers. Included with the present results are the results of Ghia, 
Ghia, and Shin [24], Schreiber and Keller [25], and Kim and Moin [6]. Listed below the flow 
quantities is the grid size used for the calculation. Good agreement among all calculations is seen 
in the lower Reynolds number cases. The discrepancies between different solutions increase at 
the higher Reynolds numbers, although the same general trend of a leveling off and then a slight 
decrease in the value of the stream-function is seen. 

To study in more detail the 10,000 Reynolds number case, the streamlines are plotted in figure 
3. The values of the stream-function contours for this plot are given in table 2. The contour levels 
plotted correspond with the values plotted by Ghia, Ghia, and Shin [24] for this case. Qualitatively, 
the plots appear to be identical. They each show secondary vortices of the same size and shape in 
the lower comers and the upper left comer. In table 3, the stream-function, vorticity, and location 
of the vortex core for all the secondary vortices for this 10,000 Reynolds number case are given 
for both the present results and the results of Ghia, Ghia, and Shin [24]. In this table, the initial 
T stands for top, B stands for bottom, R stands for right, L stands for left, and the superscript 
number corresponds to the level of the secondary vortex. Thus BR 3 refers to the third and smallest 
secondary vortex found in the bottom right comer of the cavity. Good agreement between the two 
computations is seen for this case, especially considering that the results of Ghia, Ghia, and Shin 
[24] uses over 10 times as many grid points (66,049 versus 6561). 

The convergence toward a steady-state for this problem was very good for the three lowest 
Reynolds number cases, which required less than 100 iterations and only 35 sec of computing time. 
The higher Reynolds number cases were slower to converge; the 10,000 case took 550 iterations 
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and 215 sec of computing time. The average of the computing requirements for all seven cases 
came out to 250 iterations and 100 sec of computing time. 


Flow Over a Backward Facing Step 


A second 2-D problem which has been used as a validation case is the flow over a backward- 
facing step. The challenge in modeling this problem comes from the fact that the size of the sepa- 
ration bubbles downstream of the step are very sensitive to the pressure gradient in the flow. The 
geometry used in the calculations is shown in figure 4. At the inflow boundary, a parabolic profile 
is prescribed throughout the calculation, and the static pressure is allowed to change. Two step 
heights downstream from the inflow a two-to-one expansion is encountered. The outflow bound- 
ary extends to 30 step heights downstream of the step. The ability of the flow code to predict the 
reattachment length, xl, of the primary separation bubble behind the step, as well as the separation 
and reattachment locations, x2 and x3, of the secondary separation bubble on the opposite wall, 
was tested by comparing the computed results to experimental values given by Armaly et al. [28]. 
These quantities were measured for the laminar range of Reynolds numbers, which are based on 
the average inflow velocity and twice the step height. The flow was calculated using a grid of 100 
points in the stream wise direction and 53 points in the crossflow direction. The initial conditions 
were specified to be free-stream velocity at the interior points with uniform pressure everywhere. 
For the Reynolds numbers of 100 and 200, (3 was set to 1, for the Reynolds number of 300 case, 
0.5 was used, and for the Reynolds numbers of 400 through 800, /3 was set to 0.1. The implicit 
line-relaxation process used sweeps along the primary flow direction. 

In figure 5, the quantities xl, x2, and x3 are plotted versus Reynolds number for both the 
present computed results and the experimental results of Armaly et al. [28]. Good agreement is 
seen between the two for the value of xl at the lower Reynolds numbers before the secondary sep- 
aration appears. At a Reynolds number of 400, the secondary separation bubble is present, and the 
computed primary reattachment length begins to fall off of the experimental curve. Similarly, the 
computed secondary separation points are shorter than the experimental values, although the same 
behavior is seen; that is that the secondary separation point is upstream of the primary reattach- 
ment point. The computed secondary reattachment point is seen to be close to the experimental 
values. In their experiment, Armaly et al. reported that the flow was found to be 3-D near the 
step for Reynolds numbers greater than 400, and that the 3-D effects were negligible for lower 
Reynolds numbers. These 3-D effects could explain the discrepancies between the calculations 
and experiment. 

Results similar to the present results were reported by Kim and Moin [6]. They reported a 
primary reattachment length of just under 12 step heights for a Reynolds number of 800, and the 
present result for this Reynolds number is 1 1 .48. They reported a secondary separation bubble size 
of 7.8 and 11.5 step heights for Reynolds numbers of 600 and 800, respectively. The present results 
show secondary separation bubble sizes of 7.34 and 11.07 step heights for these two Reynolds 
numbers. The similarities between these computational results give credence to the idea that the 
three-dimensionality of the flow affects the separation bubble size. 

The convergence characteristics of the code for this problem are very good. In figure 6 the 
convergence histories of the Re = 100 and 800 cases are plotted. Figure 6a plots the log of the max- 
imum residual normalized by the maximum residual at the first time step versus iteration number 
for the Re = 100 and 800 cases. Figure 6b plots the primary reattachment length xl versus iteration 
number. The Re = 100 case converges within 55 iterations and the Re = 800 case converges within 
165 iterations. The average number of iterations required for the eight different Reynolds number 
cases is 104 and the average required computing time is just under 1 1.5 sec. 


14 



Flow Over a Circular Cylinder 


As an example of an external flow problem, the flow over a 2-D circular cylinder was cal- 
culated. The grid was an algebraically generated o-grid with 100 points in the circumferential 
direction and 60 points radially. The grid points were clustered radially toward the body and the 
outer boundary was placed 10 diameters from the cylinder. The code was run and steady- state 
solutions were obtained for the Reynolds numbers of 5, 10, 20, and 40, based on the free-stream 
velocity and the cylinder diameter. The value of the artificial compressibility constant /? was set to 
50 for all the cases. At the outer boundary, where fluid was entering the domain, the velocity was 
held constant, and where fluid was leaving the domain, the static pressure was held constant. The 
line-relaxation scheme used four sweeps in both of the coordinate directions, which seemed to be 
the best trade-off of convergence versus computing time in numerical tests. 

For each case, several flow quantities of the flow were computed. Figure 7 shows a schematic 
diagram of the geometry for this flow problem along with several of these flow quantities. These 
quantities are the flow separation length measured from the rear of the cylinder in cylinder diameters 
(L se o), the angle which defines the point of separation from the body (0 sep ), the coefficient of drag 
(Cx>), the coefficient of pressure drag (Cd p ), and the coefficient of pressure at the front (C p f) and 
rear (C p r ) stagnation points. In table 4, these quantities are presented for the present calculations as 
are the numerical results of Takami and Keller [29], Dennis and Chang [30], Tuann and Olson [31], 
Braza, Chassaing, and Minh [32], and the experimental results of Coutanceau and Bouard [33] and 
of Tritton [34]. The comparisons show that there is very good agreement among nearly all of the 
results. 

The convergence of the code toward a steady- state solution for the problem of a circular cylin- 
der was found to be good. All four Reynolds number cases converged in less than 70 iterations, 
requiring an average of 21 sec of computing time. 


CONCLUSIONS 


The use of a flux-difference split upwind-differencing scheme has been applied to the artificial 
compressibility equations to solve the steady-state Navier-Stokes equations. This eliminates the 
need for any explicitly added artificial dissipation terms and any artificial dissipation coefficients. 
Although no direct comparison has been made, this code has been found to be much more robust and 
easier to run than previous applications by the authors of the artificial compressibility method which 
used central differencing plus artificial dissipation. The natural addition of dissipation through the 
use of upwind-biased stencils requires no trial-and-error adjusting of smoothing parameters as does 
the artificial dissipation. As well, it is thought that the terms on the main diagonal of the implicit 
side matrix which are not present in the central difference scheme make the scheme much more 
robust. Implicit boundary conditions based on the method of characteristics were presented. The 
accuracy of the upwind scheme has been established using three 2-D test cases. Good comparison 
was found between the current method and other methods which used many more grid points in 
calculating the flow inside a driven cavity. Results which compared well with experimental values 
were obtained in calculating the flow over a backward-facing step. The discrepancy at higher 
Reynolds numbers could be explained by the three-dimensionality of the experiment, and this was 
supported by other computational results. Finally, good comparison was found in measuring the 
flow around a circular cylinder. Perhaps the most striking feature of the current code is its ability 
to obtain steady-state solutions in a small number of iterations for most problems. Very good 
convergence rates were observed when the proper choice of the artificial compressibility constant 
j3 was made. No analytical guidelines for the choice of this parameter have been derived as yet, and 
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this is a subject of ongoing work. The extension of the current code to three dimensions is nearly 
completed and has been found to be straightforward. 
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TABLE 1 

STREAM-FUNCTION AND VORTICITY AT THE CENTER OF THE 
PRIMARY VORTEX FOR DIFFERENT REYNOLDS NUMBERS 


Re 

Present 

Ipmin (^v.c.) 

Ghia et al. [24] 

V’rmn ( ^v.c . ) 

Schreiber etal. [25] 
V’min (W„. c .) 

Kim et al. [6] 
V'm.n (u>„.c.) 

100 

-0.1030(-3.104) 

81x81 

-0.1034(-3.166) 

129x129 

-0,1033(-3.182) 

121x121 

-0.1030(-3.177) 

65x65 

400 

-0.1131(-2.296) 

81x81 

-0.1139(-2.294) 

129x129 

-0.1130(-2.281) 

141x141 

-0.1120(-2.260) 

65x65 

1000 

-0.1 171 (-2.044) 
81x81 

-0.1179(-2.050) 

129x129 

-0.1160(-2.026) 

141x141 

-0.1160(-2.026) 

97x97 

3200 

-0.1 195 (-1.904) 
81x81 

-0.1204(- 1.989) 
129x129 

- 

-0.1150( 1.901) 
97x97 

5000 

-0.1192(- 1.846) 
81x81 

-0.1190(- 1.860) 
257x257 

- 

-0.1120(1.812) 

97x97 

7500 

-0.1186( 1.846) 

81x81 

-0.1200(- 1.880) 
257x257 

- 

- 

10000 

-0.1177( 1.826) 
81x81 

-0.1197(-1.881) 

257x257 

-0.1028(- 1.622) 
180x180 

- 


TABLE 2 

VALUES FOR STREAMLINE CONTOURS IN FIGURE 3 

Contour 

number 

Value of ip 

Contour 

letter 

Value of 

1 

1.0 x 10" 5 

A 

-1.0 x 10" 5 

2 

5.0 x 10 " 5 

B 

-1.0 x 10~ 4 

3 

1.0 x 10" 4 

C 

-0.01 

4 

2.5 x 10 ~ 4 

D 

-0.03 

5 

5.0 x 10 " 4 

E 

-0.05 

6 

1.0 x 10- 3 

F 

-0.07 

7 

1.5 x 10- 3 

G 

-0.09 

8 

3.0 x 10 ~ 3 

H 

-0.1 



I 

-0.11 



J 

-0.115 
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TABLE 3 

PROPERTIES OF THE SECONDARY VORTICES FOR THE DRIVEN CAVITY AT Re = 10,000 


Vortex 

Results 

4>v.c. 

W v .c. 

X„. c . 

y v . c . 

TL 

Present 

2.418xl0 -3 

2.191 

0.0723 

0.9117 


Ghia, Ghia, and 
Shin [24] 

2.420x1 0" 3 

2.183 

0.0703 

0.9141 

BL 


1.434xl0 -3 

2.097 

0.0585 

0.1686 



1.518xl0 -3 

2.086 

0.0586 

0.1641 

BR 


3.227xl0 -3 

4.163 

0.7619 

0.0585 



3.418xl0 -3 

4.053 

0.7656 

0.0586 

BL 2 


-5.120xl0" 7 

-0.02207 

0.1416 

0.01722 



-7.757xl0- ? 

-0.02754 

0.1560 

0.01950 

BR 2 


-2.103xl0- 4 

0.3726 

0.9277 

0.07288 



-1.313xl0 -4 

0.3126 

0.9336 

0.06250 

BR 3 


-4.267xl0 -7 

-2.956xl0" 3 

0.9981 

0.008697 



-5.668xl0 -9 

- 

0.9961 

0.003900 
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TABLE 4 

FLOW QUANTITIES FOR A CIRCULAR CYLINDER 


Source 



Reynolds Number 



5 

10 


20 

40 




Lsep 


Present 

0 

0.254 

0.932 

2.29 

[29] 

0 

0.249 

0.935 

2.32 

[30] 

0 

0.252 

0.94 

2.35 

[31] 

0 

0.25 


0.9 

2.1 

[33] (exp) 

- 

0.34 


0.93 

2.13 




6 aep (degrees) 


Present 

0 

28.8 


43.1 

53.0 

[29] 

— 

29.3 


43.7 

53.6 

[30] 

0 

29.6 


43.7 

53.8 

[31] 

>6 

29.7 


44.1 

54.8 

[32] 

— 

— 


43.6 

54.5 

[33] (exp) 

— 

32.5 


44.8 

53.5 




C D (C Dp ) 


Present 

4.18(2.19) 

2.89 

(1.602) 

2.08 (1.242) 

1.549 (1.011) 

[29] 

- 

2.80 


2.01 

1.536 

[30] 

4.12 (2.20) 

2.85 

(1.600) 

2.05 (1.233) 

1.522 (0.998) 

[31] 

4.66 (2.48) 

3.18 

(1.775) 

2.25 (1.35 ) 

1.675 (1.095) 

[32] 


— 


2.18 

1.60 

[34] (exp) 

4.16 

3.06 


2.02 

1.65 




C P f(- 

-CV) 


Present 

1.847 (1.067) 

1.476 (0.755) 

1.265 (0.615) 

1.147 (0.536) 

[29] 

- 

1.474 (0.670) 

1.261 (0.537) 

1.141 (0.512) 

[30] 

1.872 (1.044) 

1.489 (0.742) 

1.269 (0.589) 

1.144 (0.50Q) 

[31] 

2.23(1.081) 

1.744 (0.773) 

1.457 (0.614) 

1.312(0.543) 
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Fig 2. Comparison between present results (solid line) and computations of Ghia et al. [17] (sym- 
bols). □ : Re=100, O- Re=400, A: Re=1000, +: Re=3200, x: Re=5000, <0: Re=7500, and 
V : Re= 10,000. (a) U-velocity component versus y on the vertical centerline, (b) V- velocity 
component versus x on the horizontal centerline. 
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Fig. 4. Geometry of backward-facing step flow problem. 



Fig. 5. Separation length versus Reynolds number for the flow over a backward-facing step. Solid 
line: computed xl, dashed line: computed x2, dotted line: computed x3, A: experimental xl, 
+: experimental x2, and x : experimental x3. 
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Fig. 6. Convergence history for the flow over a backward-facing step. A: Re = 100, x: Re = 800. 
(a) Log of the maximum residual versus iteration number, (b) Primary reattachment length 
versus iteration number. 
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